“Calhoun 


Institutional Archive of the Naval Postgraduate School 





Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 1. Thesis and Dissertation Collection, all items 


1996-09 


Use of the symmetrical number system in 
resolving undersampling aliases 


Leino, Richard E. 


Monterey, California. Naval Postgraduate School 
http://ndl.handle.net/10945/32261 


This publication is a work of the U.S. Government as defined in Title 17, United 
States Code, Section 101. Copyright protection is not available for this work in the 
United States. 


Downloaded from NPS Archive: Calhoun 


Calhoun is the Naval Postgraduate School's public access digital repository for 
| (8 D U DLEY research materials and institutional publications created by the NPS community. 
«ist sia Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS's first 


NY KNOX appointed — and published -- scholarly author. 

ia) LIBRARY Dudley Knox Library / Naval Postgraduate School 

411 Dyer Road / 1 University Circle 
Monterey, California USA 93943 





http://www.nps.edu/library 





NAVAL POSTGRADUATE SCHOOL 
Monterey, California 





THESIS 


USE OF THE SYMMETRICAL NUMBER 
SYSTEM IN RESOLVING 
UNDERSAMPLING ALIASES 


by 
Richard E. Leino 
September 1996 
Thesis Advisor: | Phillip E. Pace 





Approved for public release; distribution unlimited. 





DTIC QUALITY INSPECTED 1 














Form Approved 
OMB No. 0704-0188 





REPORT DOCUMENTATION PAGE 





Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering 

and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of 

information, sincera suggestions for reducing this burden to Washington bias aaah Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 
rina 0 Pp iect Pi 


1. AGENCY USE ONLY (Leave blank) | 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED 
September 1996 Master’s Thesis 


4. TITLE AND SUBTITLE 5, FUNDING NUMBERS 
USE OF THE SYMMETRICAL NUMBER SYSTEM IN 


RESOLVING UNDERSAMPLING ALIASES 
6. AUTHOR(S) 


Leino, Richard E. 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION 
REPORT NUMBER 





Naval Postgraduate School 
Monterey, CA 93943-5000 





9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 


11. SUPPLEMENTARY NOTES 


The views expressed in this report are those of the author and do not reflect the official policy or 
position of the Department of Defense or the United States Government. 


12a. DISTRIBUTION/AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE 


Approved for public release; distribution unlimited. 


13. ABSTRACT (Maximum 200 words) 


Two algorithms are presented which allow for the unambiguous resolution of multiple undersampled 
frequency components in a signal. Digital signal processing is usually governed by the Nyquist criterion which 
limits the amount of information that can be unambiguously stored and recovered digitally to a spectral width 
no larger than half the sampling frequency. Both algorithms resolve a spectrum beyond Nyquist by using 
additional information. The first method samples a signal more than once using a different sampling frequency 
each time. The second method utilizes a single sampling frequency which is used to sample both the signal and 
a band-limited version of the signal. When using multiple sampling frequencies, each sampling frequency 
yields a digital sequence which, in turn, has a unique spectrum when the Discrete Fourier Transform (DFT) is 
applied. The bin and amplitude information from each of the resulting undersampled spectra is then 
recombined to resolve the original spectrum. In like manner, when using a single sampling frequency the 
spectra of both the signal and its band-limited version are recombined to obtain the solution. Given a sampling 
frequency, both algorithms allow for the unambiguous resolution of a signal with a spectral width at least twice 
as large as that predicted by Nyquist. 


14. SUBJECT TERMS 15. NUMBER OF PAGES 


17. SECURITY CLASSIFICATION | 18. SECURITY CLASSIFICATION | 19. SECURITY CLASSIFICATION | 20. LIMIITATION OF 
OF REPORT OF THIS PAGE OF ABSTRACT ABSTRACT 


Unclassified Unclassified Unclassified UL 
NSN 7540-01-280-5500 1 STANDARD FORM 298 (Rev. 2-89) 





I 











Approved for public release; distribution is unlimited 


USE OF THE SYMMETRICAL NUMBER SYSTEM IN RESOLVING 
UNDERSAMPLING ALIASES 


Richard E. Leino 
Captain, United States Marine Corps 
B.S., Rensselaer Polytechnic Institute, 1990 


Submitted in partial fulfillment of the 
requirements for the degreee of 


MASTER OF SCIENCE IN ELECTRICAL ENGINEERING 


from the 


NAVAL POSTGRADUATE SCHOOL 
September 1996 


Author: 








Approved by: 
Phillip E. Pace, Advisor 


Dayid Styer =n Rea 


Herschel! H. Loomis, Jr., Cha 
Department of Electrical and Computer Engineering 







ill 





1V 














ABSTRACT 


Two algorithms are presented which allow for the unambiguous resolution of 
multiple undersampled frequency components in a signal. Digital signal processing 
is usually governed by the Nyquist criterion which limits the amount of informa- 
tion that can be unambiguously stored and recovered digitally to a spectral width 
no larger than half the sampling frequency. Both algorithms resolve a spectrum be- 
yond Nyquist by using additional information. The first method samples a signal 
more than once using a different sampling frequency each time. The second method 
utilizes a single sampling frequency which is used to sample both the signal and a 
band-limited version of the signal. When using multiple sampling frequencies, each 
sampling frequency yields a digital sequence which, in turn, has a unique spectrum 
when the Discrete Fourier Transform (DFT) is applied. The bin and amplitude infor- 
mation from each of the resulting undersampled spectra is then recombined to resolve 
the original spectrum. In like manner, when using a single sampling frequency the 
spectra of both the signal and its band-limited version are recombined to obtain the 
solution. Given a sampling frequency, both algorithms allow for the unambiguous 
resolution of a signal with a spectral width at least twice as large as that predicted 


by Nyquist. 
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I. INTRODUCTION 


A. BACKGROUND 

For a periodic discrete signal z,(n7') with a period of N samples or a discrete 
signal z(nT) that contains exactly N samples, a transform known as the discrete 
Fourier transform (DFT) is defined. The DFT transforms the N time-domain sam- 
ples to N frequency-domain coefficients where the frequencies are 27k /NT where 
k = 0,1,...,N —1 and T is the sampling period. The DFT coefficients X(k) are 
also periodic with period N. For actual computations using a digital computer, the 
continuous signal is digitized, the number of samples is limited, and the Fourier rep- 
resentation is obtained for a finite number of frequency values. The digitization of the 
signal is usually governed by the Nyquist criterion where it is assumed that the input 
signal must be bandlimited (0 < f < f;/2) before going into the analog-to-digital 
converter (ADC). The Nyquist theorem however, only places a limitation on the in- 
formation that can be derived from a single set of digitized data [3]. That is, a single 
set of digitized data limits subsequent analysis to an f,/2 bandwidth unless there is 
additional information available. With additional information, the frequency compo- 
nents f > f,/2, which appear ambiguously due to undersampling, may be resolved 


provided they are within the dynamic range of the algorithm utilized. 


B. UNDERSAMPLING 

There are several advantages to an undersampled system [4]. Among these are 
a reduction in the speed requirements on the digital section of the system, a relaxation 
on ee analog anti-aliasing filter requirements, and the possibility of extending the 
capabilities of existing systems with relatively minor redesign. Also, a good deal 


of power and money can be saved in the analog-to-digital converter section. The 





main problem in identifying the frequencies present in an undersampled ‘signal is the 
resolution of the ambiguities. 

There has been a great deal of research investigating methods to either de- 
tect aliasing or to recover undersampled signals. Baier and Furst [1] have proposed 
a, method to detect aliased frequency components by utilizing two ADC’s at slightly 
different sampling rates. The frequency components occuring in the resulting spectra 
are also slighty shifted and the difference in the two spectra allows for detection of 
the aliased frequency component. Barbarossa [2] has proposed a method which can 
recover an undersampled signal using time-frequency domain analysis and computing 
the Wigner-Ville distribution of the signal. Rader [10] has described an undersam- 
Sia technique that can recover periodic signals by reconstructing the waveform using 
a set of trial sampling periods. The trial period which yields the waveform of small- 
est variation is then considered to be the correct period and the resulting waveform 
the correct waveform. Zoltowski and Matthews [13] have devised an algorithm for 
unambiguous estimation of both frequency and phase for a single frequency utilizing 
eigenvector information. To come more in line with real-time wideband processing, 
methods based on the use of phase shift information to resolve the ambiguities in a 
single frequency undersampled signal have also been investigated [12], [11]. Finally, 
Martin and Rasmussen [6] have proposed a testbed for evaluating the impact that 
ADC imperfections would have towards limiting digital implementations of under- 


sampling. 


C. PRINCIPLE CONTRIBUTION 
Current research in the area of undersampling is limited in two principal ways. 
First, while methods exist that can unambiguously detect a single frequency, none 


of the research proposes a dynamic range to its system or algorithm for detection. 

















Secondly, no algorithm currently exists which allows for the unambiguous detection 
of a signal containing a spectrum of undersampled frequencies. 

The principal contribution of this thesis is to resolve both problems. ‘This thesis 
presents algorithms which solve for an undersampled spectrum of frequencies. The 
algorithms are based upon the mathematical properties of the Symmetrical Number 
System (SNS) [9] and the relationship between digital signals and the SNS. Each 
of the algorithms developed are simple linear solutions with dynamic ranges which 


ensure no ambiguity in the resolved spectrum. 


D. THESIS OUTLINE 

This thesis presents two algorithms which unambiguously resolve an under- 
sampled spectrum band-limited beyond Nyquist. For the first algorithm, it will first 
be shown that if a signal is sampled at two different sampling frequencies and the DFT 
is applied and normalized, two unique spectra result. Utilizing the bin and amplitude 
information of each spectrum as well as knowledge of how the real and imaginary 
frequency components alias into DFT bins, it is shown that the original frequency 
components may be unambiguously resolved utilizing a matrix multiplication. For 
the second algorithm, it will be shown that if both the signal and a band-limited 
version of the signal are processed, it is possible to resolve a spectrum of frequencies 
up to nearly twice the Nyquist rate. The band-limited version of the signal allows for 
the unambiguous resolution of frequencies below ds which is then applied to solve for 
frequencies Is <f<f,-—1 

Since the properties of the SNS are what the solution methods utimately uti- 
lize, Chapter II begins with a detailed look at the Symmetrical Number System (SNS) 
and its dynamic range for unambiguously detecting a tone frequency is introduced 


and proven. Chapter III] shows that the DFT encodes frequency information in a 








format that is exactly the same as the SNS which means that a digital sequence has 
the same properties as the SNS. Beginning to resolve undersampled frequencies, the 
properties of the SNS allow us to develop the architectures of Chapters IV and V, 
in which two and three channel architectures are introduced which unambiguously 
resolve a single frequency within the SNS dynamic range. Moving towards resolving 
multiple frequencies, Chapter VI shows that the DFT can be normalized and using 
the properties of the SNS, a method is demonstrated which resolves two undersam- 
pled frequencies in a signal. With the SNS and normalization, there is now sufficient 
information of resolve an entire spectrum of frequencies. In Chapter VII the bin loca- 
tions and normalized amplitudes are used to set up linear equations which solve for an 
undersampled spectrum with frequencies with zero-phase and no DF'T leakage. Since 
real world signals have random phase, Chapter VIII removes the phase constraint 
and utlimately presents the first algorithm which resolves frequencies with random 
phase and no DFT leakage and band-limited at twice Nyquist. Chapter IX presents 
the second algorithm which requires only one sampling frequency but has the same 
dynamice range as the first algorithm. Finally, in Chapter X, the limitations of the 
two algorithms are explored and a possible direction for resolving the limitation is 


explored. 











Il. THE SYMMETRICAL NUMBER SYSTEM 


A. DEFINING THE SNS WAVEFORM 

In the SNS preprocessing, r different periodic symmetrical waveforms are used 
with periods based on r different, pairwise relatively prime, integer lengths m1, ma, 
...,M,. For each integer, k, an r dimensional column vector A;, is formed by placing 
the value of the i*" waveform at k in the i** position, 1 < i <r. It is then desired to 
find out the largest set of vectors Ao, A1,...,A, that are distinct. This sequence of 
k + 1 vectors forms the unambiguous output of the system. This number, k + 1, is 
called the dynamic range of the system. Dynamic range refers to the fact that each 
vector up to A, is unique in the same way that frequencies up to f,/2 are unique and 
unambiguous when sampled and conventional methods are applied. 


The definition of the SNS waveform is given below. 


Definition 2.1 Let m be an integer greater than 1. For an integer h such that 
O0<h<m, define 


an = min{h,m — h} ; (1) 
This function is extended periodically with period m. That 1s, 
Qh+nm = Gh (2) 


where n € {0,+1,+2,...}, and ap, is called a symmetrical residue of h+ nm modulo 


™m. 


Let @,,, be the row vector [ao, a1,...,@m—1]. For m odd 
mm Mm 
am = oS eee aleloa eee fe 3 
- | A | | 3) 


where |x| indicates the greatest integer less than or equal to z. For m even, 
mm 
Oe NO Be ey Vee os 4 
em | 272 | 4) 


o 











Both have size 1 x m and consist of the symmetrical residues elements a;,,0 <h < m. 
This shows the form of one period of length m. From this definition it is emphasized 


that it follows for any integers, h, k, an = @n+x if and only if h = +(h + k)(mod m). 


B. DYNAMIC RANGE OF AN SNS SYSTEM 
Having introduced the SNS waveform, it is useful to introduce the dynamic 


range of an SNS system given its moduli. 


Theorem 2.1 Let m,...,m, ber pairwise relatively prime moduli, and let Ao, Ai, 


Ao,... be vectors formed by the symmetrical number system given in Definition 2. 1. 


a) If one of the moduli (m1) is even, then the dynamic range of the system 1s 
a F My J = 
M = min or I[ m+ [] m, (5) 
f=2 f=j+1 
where 7 ranges from 1 tor—1 and m,,,miz,...,™Mi, Tange over all permutations 
) es eee i 
b) If all of the moduli are odd, then the dynamic range of the system 1s 
: na Tae 
M = min 5 [] ™ + 5 [I] me? (6) 
l=] f=j+1 
where j ranges from 1 tor—1 and mj,,Mig,...,mi, Tange over all permutations 


OF Ae. oe 


1. Example of the SNS Dynamic Range 


Let m, = 4, m2 = 3, and m3 = 5. We must minimize the set of values: 


{a + mama, — mg + m3, img + ma} = {17, 11, 13}. 


The dynamic range is the minimum value of this set, 11, as is verified in 


Table 2.1. We see that A1; = Aj, and Aj, is the first repetitive vector. In other 


6 














| (0 12 3 45 6 7 8 9 10 11 12 13 --. 





Table 2.1: Integer Values for the m; = 4, mz = 3, and m3 = 9 


words, Ap,...,Aio is a set of 11 distinct vectors and, Ao,...,Ai, is not a set of 
distinct vectors. This shows directly that the dynamic range of this system is 11. A 
repeated vector is called an ambiguity of the system. In the above example Aj, is the 


first ambiguity of the system. 


C. PROOF OF THE SNS DYNAMIC RANGE 
Before starting the proof of Theorem 2.1 it is necessary to express the general 
idea that will be used. Suppose that mj ,...,m, are r positive pairwise relatively 


prime integers. Also suppose that 


1 1 
Oh “h+k 
ap Chk 
Ay = . = ‘ = Antk | (7) 
T T 
Qn Qhik 


so that there is an ambiguity at the position h+k (since the vector is indistinguishable 
from that at position h), where h > 0 and k > 1. 

The key to the proof is that fact that A, = Anis, if an only ifh = t(h + 
k) (mod m,) for 1 <i <r. Thus the proof involves systems of linear congruences and 
the Chinese Remainder Theorem [7]. The goal is to find the least value of h + k that 
is an ambiguity. In general all permutations of subscripts 1, 2,...,r must be worked, 


but for the sake of a clearer presentation the subscripts are not permuted. 








1. Proof of Theorem 2.1a. 


Let m , = 2m be even. 


Case I: 
h=(h+k)(modm), Il<i<j (<7) (8) 
h.=—(h+k) (mod m,), j+l<i<r (9) 

Case II: 
h =-—(h+k) (mod m,), L<i<j (10) 
h=(h+k) (mod m,), g4+1<i<r (11) 


Note that the two cases are opposite in terms of the plus or minus signs in the 
congruences. These would not be separate cases except for the fact that m, is the 
even modulus. 

a. Case [ 

Suppose congruences (8) and (9) are true, i.e., k = 0 (mod m,) for 
1<i< j. Since the moduli are relatively prime in pairs, k = 0 (mod Te, m,;). Note 
that k is even since m, is even. It follows that k = a][{_, m; for some a € {1,2,...} 


SO 


k j 


Continuing on for the moduli m,, the condition i > 7 + 1 is examined next. From (9) 
k 
=—5 (mod m;) (13) 


for 7 +1<i<-r. Since these m; are pairwise relatively prime, the Chinese remainder 
theorem guarantees that there is a unique solution h (mod []j_,,,m:) to this system. 
Therefore, there will be exactly []/_, m; solutions (mod []f_, m;). Let b’ be the least 


integer such that 


: k 
b’ II mMm-->0. (14) 
i=j+1 2 | 











Note that if we set 


3 k; 
hy = 0 I] Mi =, (15) 
i=jtl 2 
then 
0<hi< [J m (16) 
i=j+l 
and 
k 
h,= a (mod m;) (17) 


for j+1<i<vpr. This h, is thesolution (mod lis 4,™;). All solutions (mod J[]j_; ™m;) 
are of the form 


s 
h=hi+b [] m (18) 
i=j+l 


for 0 < 6 < []2., m;. Using (15) the ambiguities, h +k, will have the form 


3 Tr Tr - Tr k 
h+k=hit+b [J mt+k=0' [[ mt+d I] mt+- (19) 
i=j+l i=j+1 i=j+l 2 
Using (12) 
j r 
h+k=am|[m +b {] Mm; (20) 
i=2 i=j+1 


where b= b'+ 5b {1,2,...} anda € {1,2,...}. 

b. Case II 

Although the particulars of this case are different from those of Case I, 
the line of argument is the same, and the result is exactly the same: 

j r 
h+k=am|[m+6 [J m (21) 
i=2 i=jtl 

where a and b are positive integers. 

Now the two cases have been treated and it is evident that all the am- 
biguous values h +k, have the form (20) or (21) and these forms are indistinguishable. 


The least number of this form (or the first of many ambiguities) is found by setting 





a=b=1. Next it is shown that 
j r 
i=2 i=j+1 
is one of the ambiguities. 


Suppose 
3 i: 
i=2 i=j+l 


We revert to Case I and set a = 1 in (12) so 


j 
=m [m, 
/ ie, 


Then, in (14) 0b’ = 1, so 
k 
Ay = I M4 — 9 
1=j+1 
and 


| r 
i=2 i=j+l 





(22) 


(23) 


(24) 


(25) 


(26) 


is an ambiguity. That is, the vector corresponding to h, is the same as the vector 


corresponding to hi +k. 
The alternative is that 
r j 
Il mi,<m™ I] M4. 
t=7+1 1=2 


Now we use Case II with a = b= 1. We can set 


k ¥ i 
=, <== I] M5 
2 t=j9+1 
‘Then 
j k 
hp=m|[m-— 
1 Mm Il M7 9 


is positive, and 


j r 
1=2 ix=jtl 


10 


(27) 


(28) 











is an ambiguity. 


Recall that 


1 1 
a8 anne 
ap, Qhik 
= (31) 
i T 
ap Ahik 


if and only if h = +(h +k) (mod m,), i = 1,...,r. There are 2” possibilities for 
choice of plus or minus signs in these congruences. For 7 fixed we have looked at two 
possibilities and found the minimum ambiguity. All others can be found by varying 


7 and by permuting the order of the moduli. Thus 


ef) r 
M=min(m|[[m,+ [] mi, (32) 
l=2 l=j+l 
is the minimum ambiguity. 


2: Proof of Theorem 2.1b. 


Let m,...,m, be r odd pairwise relatively prime natural numbers. Recall 
that : : 
“h “htk 
Gh | | Shek | (33) 
ah, Dnt k 


if and only ifh = +(h+k) (mod m;) for1 <i <r. Suppose h = (h+k) (mod m,) for 
L<i<jandh=-—(h+k) (mod m;) forj7+1<i<r. The first set of congruences 


implies that k = 0 (mod m;),1<i< j. Thus 





j 
k=al[m (34) 
1=1] 
for a € {1,2,...}. The second set of congruences implies that 2h = —k (mod m,) i.e., 
Mm; 1 


forj7 +1<1i< r. Again the proof is split into two cases, this time depending on 


whether k is even or odd. 
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a. Case [ 


Suppose & even. Then 
k 
h= =e (mod m,) (36) 
forj7 +1<121<pr. Let db’ be the least integer such that 


: k 
hy=o' [| m-=2>0 (37) 
i=j+1 2 


Note b’ € {1,2...}. This is the unique solution (mod []j_;.,; ™m) to the set of 


congruences (36). Every solution to the original problem is of the form 


T ms Tr k 
h=b' [J m+ [[ m-= (38) 
i=j+l | i=j+1 2 
where 
» @ 
0<b<|[m,. (39) 
i=1 


By (34) and (38), all ambiguities A +k that arise from these values of h are of the 


form 


j Pe oe 
h+k=-|[m+- [] (40) 
2 jy 4 init] 


where b = 2(b' + b) is an even positive integer. Note that k is even if and only if a is 
even. 

b. Case Il 

Suppose that k is odd. As in Case I, all ambiguities h + k that arise 
from these values of h are of the form 

h+k=S][m+5 TT mi (41) 
i=l i=j+l 

where a and b are odd positive integers. 

Combining Cases I and II, note that the smallest value of (40), (41) is 


found when a = b = 1. Furthermore, when a = b = 1, and Tv, mi < Tiaj41 7 this 
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is an ambiguity: Let 


ee 1 


i=j+l i=l 
and k = [[Z_, mj. It is easily checked that h = h+k(mod m,) for 1 <i <j and 
h = —(h+k) (mod m,) for j+1<i<r. When J]j_.,, mi < [I-77 we must solve 
the permuted problem: h = h+k(mod m,) for j+1<i<randh = —(h+k)(mod m,) 
forl <i <j. This time k = []f,,,mi, and h = (1/2) Ty mi — (1/2) Thay mi 
give the ambiguity at h +k. Since the solution to the problem involves looking at 
all permutations this gives the minimum ambiguity. This completes the proof of 
Theorem 2.1b. 

Theorem 2.1 provides the basis for predicting the dynamic range of 
any system which encodes information in a format identical to the SNS. In fact, in 
the next chapter it is shown that the process of digital sampling encodes frequency 
information into the SNS format which allows the DFT to be treated identically to 


the SNS. 
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III. DFT AND SNS RELATIONSHIP 


A. DIGITAL SAMPLING AND ALIASING 

Consider a single frequency signal sampled at two different sampling frequen- 
cies. Digital uniform sampling of an analog waveform with frequency f produces a 
discrete waveform which is symmetrical about the sampling frequency, f,/2. Assume 
for this system that the two sampling frequencies are fs; = 10 and fsg = 11. After 
sampling, an analog input signal, z(t), becomes a discrete sequence, z(nT). This 
periodic sequence has a digital frequency given by w = 2a(f/f,). A signal with dig- 
ital frequency 0 < w < 7 is indistinguishable from a signal with digital frequency 


nn <w<(n+1)r,n=1,2,3,..., an effect known as aliasing. 


B. DIGITAL FREQUENCY MAPPING 

The digital frequency of a sampled sinusoid can be mapped into the z-domain 
as shown in Figure 3.la. For simplicity assume a sinusoid x(t) = 2cos27 ft and after 
sampling 


a(n) = 2coswn = e)”" + eI" . (43) 


If f = f,/4, this corresponds to w = 1/2. If f = f;/2, w = m. Since the signal is 
real, the signal appears in complex conjugate pairs on the z-plane. For frequencies 
between f,/2 and f,, the frequencies map back to their conjugate on the upper half 
of the complex plane. If the frequency is increased beyond f,, a full trip is made 
around the unit circle and the mapping repeats. Figure 3.1b illustrates the mapping 
with each triangle representing a full rotation around the unit circle in the z-plane. 
The abscissa represents the input analog frequency while the ordinate represents the 
digital frequency mapping. Note that an infinite number of analog frequencies will 


map into each digital frequency 0 < w <7. 
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Input Analog Frequency 
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Figure 3.1: a)Z-plane mapping of an input analog signal. b) sampled 
frequency output. 
C. DFT BINS AND THE SNS 
Recall that the DFT is given by: 
1 * 
X(k) = So a(nje FOF") — ke =0,...,N-1. (44) 
n=0 
Application of the DFT to x(n) yields a discrete spectrum where X (k) is the energy 
contained in the signal at each digital frequency w = 27k/N. The discrete spectrum 


X(k) has N indices with the digital frequency of each index given by: 


0, ee a calc] 2 (N/2 +1) | oN = 4) oa) 


is yttt, ; for N even , 
N N N N N 


(45) 
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1 N —1)/2 N+] 
0, De eee id 2 IE a Oe 


(N — 2) 


(N ~ 1) 
N N 


for N odd. 





2 

(46) 
The analog frequency corresponding to each index is obtained by multiplying each 
value by f,. Since signals with digital frequencies in the range 7 < w < 27 are indis- 
tinguishable from signals with digital frequencies 0 < w < 7, the digital frequency of 


each index can also be written as: 


1 N/2 N/2-1 2-4. SD 
0, 27—,---, ccd a cate --+, Qr—, 2a — for N even , (47) 
N N N N N 
and 
1 N/2 N/2 Z 1 
0 an ree, an NUE at ree, one an for N odd , (48) 


where |x| represents the greatest integer less than or equal to x. More simply, the 
spectrum X(k) resolves into N integer indices and incoming signals will map into 


unique bins: 


NN 
ote Sab ed for N even , (49) 
and 
NN 
ot. LEI be 2d for N odd . (50) 


That is, since frequency indices greater than N/2 are redundant for real signals, 
the highest unaliased frequency that can be observed corresponds to the N/2 index. 
From the above discussion, it is clear that the DFT maps real signals naturally into 
the symmetrical number system [9]. In this case, the modulus described in [9] and 
in Chapter II is the sampling frequency, f,. The number of indices N is given by 


N = f,T;, where Tz is the total sampling time. 


D. ILLUSTRATING THE RELATIONSHIP 


Figure 3.2 illustrates the DFT mapping for two channels where f,; = 10 and 
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feo = 11 for input frequencies f = 0 to 23. In this case, T, = 1 so N; = 10 and 


No = 11. In the figure, the abscissa corresponds to the incoming frequency while 


the ordinate corresponds to the bin the signal is resolved into. Table 3.1 displays the 
input frequency and the resulting DFT bin for each sampling frequency. Note that 
frequencies resolve as described in (49) and (50). By considering both channels, it is 
possible to unambiguously resolve signal frequencies in the dynamic range determined 
by the SNS, (0 < f < 15). 

Simply knowing that the DFT is in the SNS format and what the dynamic 
range is for any given sampling frequencies would be useless unless the information 
can be used to resolve undersampled frequencies within the dynamic range. The next 
two chapters detail two and three channel architectures as well as solution methods 
which will resolve a single undersampled frequency band-limited within the dynamic 


range of the SNS system. 
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Figure 3.2: a)DFT mapping for input frequencies f=0 to 23 for a) fs; = 10 
and b) ae = 11. 
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Table 3.1: Input Frequency and Resulting DFT Bins for 2 Channel Exam- 
ple 
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IV. RESOLVING THE TWO CHANNEL CASE 


A. THE GENERAL SETUP 

Figure 4.1 shows the block diagram of a two-channel receiver architecture to 
determine a single frequency f. In this architecture the ADC sampling frequencies 
fs; and f. are relatively prime. The DFT outputs are thresholded to detect the 
frequency bins. The frequency bins a; and ag are then used by the SNS-to-decimal 
algorithm to determine the frequency of the input signal. Let m; = fs; and m2 = fs2 
and suppose that the incoming frequency, f (unknown) lies within the dynamic range 
M of the SNS system (5), (6). Thus, f = +a;(mod mj) and f = +a9(mod mz). For 
each of these congruences either the plus or the minus is correct, but it is not known 
which. Thus, there are four sets of two equations: 


i) 


a,(mod m4), 
ao(mod ma), 


ay (mod m1), 
—aj(mod mo), 


ii) 


It All 


—a,(mod m1), 
—a (mod m2), 


iii) 


—a ,(mod m3), 
ao(mod m2) . 


iv) 


i a a as 
It UH 


in 


The Chinese remainder theorem [7] guarantees that each of these has a unique solution 
modulo m ;m2, and Theorem 2.1 guarantees that exactly one of these solutions lies 
within the dynamic range of the system and this is the value of f. In fact, it is only 
necessary to solve i) and ii), at most, because the solutions to iii) and iv) are the 
negatives of the solutions to i) and ii), respectively. 

Recall that in the standard statement of the Chinese remainder theorem we 


wish to solve for f, where f = a;(mod m;), 1 < i < r, and the m, are pairwise 
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SNS-to-Decimal Algorithm 


Fun ction 





Figure 4.1: Block diagram of a two channel receiver architecture to deter- 
mine a single frequency / 


relatively prime. The theorem states that there is a unique solution modulo M = 


mym2--:m,. A standard method of solution is to find integers b; such that 


Mb 





_ 1(mod m;) , (51) 
74 
2 =1,2,---,r, in which case the solution 
f =< M bya;/m, a M boaq/me ere eae Mb,a,/m,(mod M) : (52) 


Returning to the two channel case, note that the values of b; and bg depend 
only on m, and mg, and not at all on +a, or tag. Thus it may be assumed that 
the constants cy = m2b1(= Mb,/m 1) and cg = myb2 are known, and that the SNb- 
to-decimal algorithm only needs to evaluate +c ,a1 + coa2 modulo M(= m,mgq), and 


pick the one value that lies within the dynamic range. 
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B. SET-UP FOR A PARTICULAR TWO CHANNEL RECEIVER 


Let m, = 330 and m2 = 337. Then 
b} = —47 solves 337b, = 1 (mod 330) , (53) 


and 


bo = 48 solves 330b2 = 1 (mod 337) . (54) 


Thus, cy) = mob, = —15, 839, and co = mjibo = 15, 840. Also, M = 111, 210, and, by 
Theorem 2.1, the dynamic range of the system is m1/2+m 2 = 502. These coefficients 


are hardwired and do not need to be recomputed. 


C. SOLVING A PARTICULAR CASE - THE WORK OF THE ALGO- 
RITHM 
Suppose that the frequency of the incoming signal is f, 0 < f < 501, and 
a; = 59 and ag = 66. By Theorem 2.1, f is the unique solution to +(—15, 839)(59) + 
(15, 840)(66) modulo 111, 210 that lies in the interval [0,501]. First try (—15, 839) - 
59 + 15,840 - 66 = 110, 939; this is congruent to —271 modulo 111,210, so case iii) 
solves f = 271. 


It is emphasized that this is the only computation the algorithm needs to 


compute. 


D. SPECIAL CASES 


In special cases there are even quicker solutions than given in section C. For 


‘example, suppose m, = 2p and m2 = 2p +1. By Theorem 2.1 the dynamic range of 


the system is m,/2 + m2 = 3p +1. Frequencies within the dynamic range will fall 


into bins as follows: 


p+i terms p terms . p terms 
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The top bin represents the sampling frequency of m; = 2p and the bottom m2 = 2p+l. 
Now suppose the incoming frequency, f, is resolved into bins a; and ag, respectively. 


From the above it is clear that: 


a1 if Qa, = a2 
f = mm 1— a, if a, = a2—- 1 (56) 
my 1+ ay if ay =agtl1 


This makes the solution much more straightforward. 
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Figure 5.1: Block diagram of a three channel receiver architecture to de- 
termine a single frequency f 


Or 


f = +2524 70+120 (mod 210) , (58) 
but 252 (mod 210) = 42 (mod 210) so: 
f =+42470+120 (mod 210). (59) 
First, 
f = 42+ 70+ 120 = 232 = 22(mod 210) , (60) 


a value that can be discarded, because it and its negative are out of the dynamic 
range. 


Second, 


f =42+70 — 120 = —8(mod 210) . (61) 
Although —8 is out of range, the negative, f = 8, is in the dynamic range, so that 


f = 8 is the correct frequency value. 
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VI. RESOLUTION OF TWO UNDERSAMPLED 
FREQUENCIES 


A. DFT AMPLITUDE NORMALIZATION 
Consider a unit amplitude signal containing a single integer frequency. After 
sampling and forming an N point DFT X(k), the k™ bin corresponding to the digital 


frequency of the signal will have an amplitude A(k) given by: 


N eS 0. Or Vans 2 


ld a N/2 elsewhere (62) 


Dividing X(k) by A(k) will produce a normalized spectrum. 

Figure 6.1 shows two DFTs of an undersampled signal containing two fre- 
quency components with unknown amplitude and frequency. Figure 6.1a displays the 
frequency spectrum for f,; = 10 while Figure 6.1b corresponds to fsg = 11. The two 
frequencies are known to lie within the dynamic range of the SNS system which in 
this case corresponds to f < 15. It is difficult to determine which two SNS pairs are 
valid by looking at the spectra of Figure 6.1. That is, either the bin pairs A are 


correct or 3 : A are correct. 


B. BENEFITS OF NORMALIZATION 

The benefits of normalizing the DFT are shown in Figures 6.2a and 6.2b where 
each of the original spectra have been divided by Ai(k) and Ag(k) respectively. In 
this case N; = 10 and No = 11. It is now easy to see which magnitudes “match-up” 
and that 3 A is the correct set. Using the special case of two SNS moduli presented 
in Chapter IV, m; = 2p and mg = 2p + 1 where p = 5, it is determined that AER. 


corresponds to f=9 and f=14 with amplitudes 0.7514 and 0.6515 respectively. 


Theorem 6.1 Let f.; and fso be any two relatively prime sampling frequencies. Con- 


sidering a signal containing two unknown integer frequencies with distinct amplitudes, 


a 








Bin 


b) 
Figure 6.1: DFT for two unknown frequencies for a) fs; = 10 and b) fs = 11. 


the N-point normalized DFTs reveal two unique bin pairs which, in turn, reveal the 


two frequencies provided they are within the dynamic range of the SNS system. 


C. SPECIAL CASE 

The solution has a special case when the two frequencies alias to the same bin 
for one of the sampling frequencies. For example, incoming frequencies f = 3 and 
f = 13 would both alias into bin 3 for f,; = 10 but into bins 3 and 2 for f,2 = 11. 
Normalization of the DFT will still produce accurate results as the magnitudes at 


bins 2 and 3 for f,2 = 11 will add up to the magnitude of bin 2 for f,; = 10 and the 
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Figure 6.2: Normalized DFT for two unknown frequencies for a) f;; = 10 
and b) f.2 = 11. 


bin pairs H A will yield the correct frequencies. When two frequencies are received 
with identical amplitudes, Theorem 6.1 does not apply. 

Having shown that normalization allows us to “see” which bin pairs “match”, 
normalization can also be used to measure how much energy went into each bin. If 
two frequencies land in the same bin, the normalized energy should essentially add 


together into that bin. This information along with bin location is used in the next 


chapter to solve for an entire spectrum of frequencies. 
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VIL. INTRODUCTION TO RESOLVING MULTIPLE 
FREQUENCIES 


A. SETTING UP THE SOLUTION 

The normalization of the DFT amplitude introduced in the previous chapter 
leads to the consideration of a method to resolve any number of zero-phase frequencies 
which have no DFT leakage. By zero-phase frequencies with no leakage it is meant 
that the DFT values will be all real and that each resolves its entire energy into one 
DFT bin with no leakage. 

Consider a bandlimited signal with frequency components that can range from 
0 to 10 Hz. The amplitude and frequency of any component is unknown except that 
they are integer frequencies. Let xo,21,---,219 represent the amplitudes of each 
frequency component from f=0 Hz to f=10 Hz. Assume that the signal is sampled 
at f.4=10 Hz and f,2=11 Hz and that the DFT is applied and normalized. ‘T'wo 
normalized spectra will be formed, each containing undersampling ambiguities for 
f > 5. Let the notation bj represent the value of the normalized DFT coefficient 


for bin i, f; =m. From knowledge of aliasing it is clear that: 


bii0 = Zi tX 
bo10 = Xot+Ze 
b310 = %3+27 
b4i0 = 4+ 2% 
b510 = 25 
b5.11 = LS + L6 (63) 
bay. = X4 +27 
6311 = 23+ 28 
bo011 = XQa+2Q 
b111 = 21+ LXj0 
boi. = Zo 


1 and f = 9 alias into bin 1, while f = 2 


In other words, for f,; = 10, f 


and f = 8 alias into bin 2, and so on. This forms a set of linear equations that can 
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be written as: 


== Oood’ddadoo 2 © © 
COrrocdcn”nrdononrcdncennoo 


Oo rr COCO COCO CO OF & 


amo orco°gnccco cr © & 


oo oorcod OF CO OC O&O 


Mo oo cdrrooc co oO 


Oo ooor OF OC CO © 


oo oorcaoo or & © 


QO 1 
1 0 
0 0 
0 0 
0 0 
0 0 
0 0 
1 0 
0 1 
0 0 
0 0 


Om Cc © @ © © © © © © 


61,10 
6210 
63,10 
64,10 
b5 10 
O5 11 
643 
63.11 
bo 11 
6111 
60,11 





(64) 


This is of the form Az = b where zx represents the actual amplitudes of each frequency 


component and b represents the total amplitude that is aliased into each bin for the 


two sampling frequencies. From (63) it is easily seen that the values zp through z19 


may be uniquely solved for. To begin with rp = bo1; and z5 = b5.19. Then the values 


Zo, La, L7, etc., may be solved for successively. 


LO 
Ly 
L2 
L3 
L4 


L5, = 


L6 
L7 
£8 
Lg 
L10 


oOo OO OO & © 


© 


—] 
—1 
—] 


0 
=] 
—] 
—] 
—] 

0 


— pe pt 


Thus, A is nonsingular and z = A7~‘b: 


mHSemaord°’rcidcicddre © 


Ri © OM ©: CoO: ©. Dm 


3 SG 6 OO "1 


D110 
b2 10 
63,10 
6410 
bs 10 
b5,11 
b4 1 
63.11 
b211 
bi 
boi 


(65) 


From (65), it is now possible to reconstruct the original spectrum based upon 


the normalized bin values for the two undersampled spectra. 


Where before only 


frequencies from 0 to 5 Hz could be recovered, 65 allows integer frequencies from 0 to 


10 Hz to be unambiguously recovered. It can be noted that the A matrix is 11 x 11 


and that the additional equation for bin bp9i9 was not included. It turns out that 


including this bin does not add any additional information. 
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B. EXAMPLE 


Suppose we have an incoming signal, 


s(t) = 29+ 21cos(27t) + xq cos(4rt) + x3 cos(67t) + x4 cos(87t) + 
ts cos(10t) + xg cos(12mt) + x7 cos(14rt) + zg cos(167t) + 


tg cos(18rt) + x19 cos(20zt) (66) 


All z, are unknown. The signal is sampled at f,; = 10 Hz and fs. = 11 Hz 
and the DFT is applied to both with N,=10 and No=11 points respectively. This 
produces two spectra where any frequencies from 6 to 10 Hz are aliased into bins from 
0 to 5 Hz. Both DFT’s are normalized and the results for bins 0 through 5 are shown 


in Figure 7.1a for f,;=10 and Figure 7.1b for f,g=11. 


From the values in each bin it is possible to form the b vector 
b* = [0.894 1.211 2.072 1.960 0.384 0.903 2.762 1.283 1.526 0.878 0.218] (67) 
and set up the equation z = A~‘b which yields: 
gt = [0.218 0.047 0.679 0.751 1.441 0.384 0.519 1.321 0.532 0.847 0.831] (68) 


The solution (68) can be compared with the actual values for the x, which is done in 
Figure 7.2. 
Figure 7.2a shows the normalized DFT spectrum that results from sampling 
the same signal at 20 Hz. Figure 7.2b plots the solution (68). The solution is exact. 
The solution method presented so far does little good if only signals with zero- 
phase can be resolved. In the next chapter, it is investigated how the imaginary 
components of a signal alias into bins and the results are used to generalize the 


solution for signals with random phase. 
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Figure 7.1: Normalized DFT for frequencies from 0 to 10 for a) fs; = 10 
and b) fz = 11. 
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Figure 7.2: Frequency components and amplitudes: a) actual and b) cal- 


culated 
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VIII. RESOLVING MULTIPLE FREQUENCIES: FIRST 
ALGORITHM 


A. CONSIDERING RANDOM PHASE SIGNALS 

Thus far attention has been restricted to signals with zero phase. Conse- 
quently, the DFT yields an all real spectrum. When the frequency components have 
random phase, however, a modification to (65) must be introduced. 

1. Illustration of Random Phase by Example 

Consider a signal with integer frequency components from 0 to 10 Hz given 


by: 
s(t) = 29 + x1 cos(27t + $1) + v2 cos(4at + G2) +--+ + 219 Cos(207t + G10) — (69) 


where ¢,, represents some random phase for each frequency component. If this signal 


was sampled at 20 Hz, DFT applied and normalized, bins 0 through 10 would have 


the values: 
Co ZO 
Cy x1 cos(¢1) + 7x1 sin(¢) 
c2 | = | 22C0s(d2) + jx2sin(¢2) (70) 
C10 £19 COS($10) + 7X10 Sin(¢10) 


where c; represents the complex value of bin n and c, = a, + 7b. If the same signal 
was sampled at 10 Hz, the DFT applied and normalized, the real part of each bin 


will be given by: 


21,10 r1 cos(¢1) + 29 Cos(¢g) 

2210 55 cos(¢2) + Ze cos(¢g) 

a3io | = | £3c0S(¢3) + £7 cos(¢7) (71) 
2410 x4 cos(o4) + x6 cos( de) 

25,10 5 c0s(¢s) 
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2. The Imaginary Component 


Consider the imaginary component of each bin. If the DFT is treated as a unit 
circle as in Figure 3.1, then for f,=10 Hz, all frequencies below 5 Hz lie in the top 
half of the unit circle and appear, without aliasing, in bins 0 through 5 of the DFT. 
All frequencies from 6 to 10 Hz, however, lie in the bottom half of the unit circle and 
alias back to the top half into bins 0 through 5 of the DFT. Since the real part of the 
DFT is an even function, the real components of the signals add together into bins 0 
through 5. The imaginary component of the DFT is an odd function, however, and 
while the imaginary component of frequencies from 0 through 5 Hz add into bins 0 
through 5, the imaginary components of frequencies from 6 through 10 Hz subtract. 


In other words: 


b1 10 21 sin(di) — rg sin(dg) 

bo 10 xq sin(¢2) — rg sin(¢g) 

b3i0 | = | x3sin(¢3) — x7 sin(¢7) (72) 
ba 10 x4 sin(d4) — x6 sin( 6) 

bs 10 5 sin(¢s) 


B. SOLVING FOR REAL AND IMAGINARY COMPONENTS 
Based on how the imaginary components of the signals alias into the bins, it 


is again possible to form a set of linear equations. 


6110 Z1 — Zo 

62.10 22 — Z8 

b3 10 23 27 

6410 Z4— 26 

6510 45 

b5.11 = Imag Z5 = ZE (73) 
6411 Z4—- 27 | 
6311 23 —- Zg 

6211 Z2— Zg 

b111 Z1 — 210 

bo11 Zo 
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or 


010000 0 0 0-1 OF % bd 
001000 0 0-1 0 0 Y b2 10 
000100 0-1 0 0. 0 Yo b3 10 
000010 -1 0 0 0 0 Y3 6410 
000001 0 0 0 0 0 Y4 bs 10 
000001 -1 0 0 0 9Q Ys | =] O51 (74) 
000010 0-1 0 0 0 Y6 bai 
000100 0 0-1 0 90 Y7 6311 
001000 0 0 0-1 0 Y2 bo 14 
010000 0 0 0 90 =-!1 Yo bi 
100000 0 0 0 0 90 Yi0 60,11 

where Z, = 2n(cos(@n) + 7Sin(gn)) = Xn + jY, and may be thought of as the 


normalized real and imaginary components of the signal for frequencies 0 through 10 
Hz. The solution (74) may be written as A;Y = b. A; has an inverse A; ‘ which can 
be seen in the same manner as A had an inverse in Chapter VII and (74) may be 


written as Y = Az'b or: 


Y 00000 0 0 0 0 0 17 Ff bi 
Y; 11111 -1 -1 -1 -1 OO] | dato 
Yo 01111 -1 -1 -1 O 0 0} | d3i6 
Y3 00111-1 -1 0 O O OF | bao 
Y4 00011-1 0 0 0 O 0} } b5,10 
¥y/=/00001 0 0 0 0 0 Of | b51 (75) 
Ye 00001-1 0 0 0 0 Of | bai 
Y; 00011 -1 -1 0 0 O O0f |] bdsn 
Yg 0011211 -1 -1 -1 O 0 Of | bon 
Yo O 1111 -1 -1 -1 -1 +O Of] fin 
Yi0 121212112 #21 -1 -1 -1 -1 «04 | bon 
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which solves for the imaginary component while the real component is solved as in 


(65): 
Xo 0 0 0 0 0 0 0 0 001 21,10 
X41 1 1 1 l 1 —-1 -1 -1 -1 0 0 Q2.10 
X9 0 1 1 1 1 -1l -1 -l 0 0 0 23,10 
X3 0 0 1 1 1 -1 -l 0 0 0 0 &4,10 
X4 0 0 0 1 1 -1 0 0 0 0 0 25,10 
AS = 0 0 0 0 1 0 0 0 0 0 0 Q5,11 (76) 
X6 0 0 0 0 -l 1 0 0 0 0 0 G4 11 
ve 0 0 0-1 -1 1 £212 40 40004] 4I agus 
X8 0 0 -l1 -l1 -l 1 1 1 0 0 Q 22,11 
X9g 0 -1 -1 -1 -l 1 1 1 1 0 0 Q1,11 
X10 —l] -l | —-1 -l -l 1 1 1 1 1 0 Q0,11 


In summary, when the incoming signal has components with random phase, both the 


real and imaginary components may be solved separately by the equations above. 


C. THE PIVOTAL FREQUENCY 


An inspection of the specific case presented illustrates that the bin information — 


at f = 5 Hz can be thought of as pivotal in that its accuracy affects the solution for 
the entire spectrum. Unfortunately, while the phase for the frequency component at 
f = 5 Hz (in this case) is random, the phase of Zs is not. In all cases for an N-point 
DFT, where N is even, the imaginary component of the N/2 bin will be zero. Thus, 
although the real part of the solution will be accurate, the imaginary component 
will not. In order to solve this problem with the DFT for the simple case presented 
(fo: = 10 Hz and f,2 = 11 Hz), the signal is band-limited from 0 to 5 Hz, sample at 
11 Hz, apply the DFT and normalize. The value of this DFT at bin 5 will be the 


correct value which can replace the values in (75) and (76) for 6519 and as 10. 


D. EXAMPLE 


Suppose there is an incoming signal, 


s(t) = 0.218 + 0.047 cos(2mt + 5.1369) + 0.679 cos(4mt + 4.7491) + 
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0.751 cos(6rt + 2.9044) + 1.441 cos(8xt + 5.9776) + 
0.384 cos(107¢ + 3.9756) + 0.519 cos(12mt + 2.7604) + 
1.321 cos(14rt + 5.1817) + 0.532 cos(16mt + 4.3290) + 


0.847 cos(18rt + 4.4121) + 0.831 cos(207t + 6.2024) (77) 


Assuming these amplitudes and phases are unknown, the signal is sampled 
at fs) = 10 Hz and f.so = 11 Hz, the DFT is applied and normalized. Figure 8.1 
displays the real and imaginary components of the normalized DFT for f,; = 10 Hz. 


Figure 8.2 displays the real and imaginary components of the normalized DFT for 


Fs2 = 11 Hz. 


Real Grid) 


Imaginary(b, 19) 





Figure 8.1: Normalized DFT components at f,; = 10 Hz: a) real compo- 
nent and b) imaginary component. 
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Figure 8.2: Normalized DFT components at f,. = 11 Hz: a) real compo- 
nent and b) imaginary component. 


Note from Figure 8.1b that bs 39; = 0 which is the imaginary component of 
the fifth bin for f,;=10. This is not the actual value and the actual value must be 
obtained by bandlimiting the signal between 0 and 5 Hz, sampling again at 11 Hz, 
and applying/normalizing the DFT. This yields the correct value for 6519 = —0.2844. 


Now we may form the two vectors: 


a’ = [-0.2312 —0.1741 —0.1325 — 0.8925 -— 0.2580 -— 0.7398 1.9717 


0.9290 — 0.2256 0.8476 0.2180] 
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and 


by = (0.7663 — 0.1852 1.3546 -—0.6266 -— 0.2844 -— 0.4775 0.7446 


0.6698 0.1306 0.0242 0.0000] : (78) 


Solving X = Avta and Y = Aj’b yields: 


Zz. 0.2180 
Zi 0.0194 — 70.0428 
Z, 0.0249 — 70.6785 
J, —0.7300 + 70.1765 
Zi 1.3742 — 70.4335 
Zs | = | —0.2580 — 70.2844 (79) 
LAr —0.4817 + 70.1931 
Wie 0.5975 — 71.1782 
oe —0.1990 — 70.4934 
Zs —0.2505 — 70.8091 
ae 0.8283 — 0.0671 


A simple check reveals that this is the correct result. 

This simple example helps to illustrate the solution method. Systems have 
been simulated in which fs; = 1000 and f,. = 1001 and exact results are attained in 
every case. One of the drawbacks of the solution is the need to band-limit the signal 
and sample a third time at f.2 in order to recover the phase of the frequency compo- 
nent at fs,/2. An alternative method might be to simply not transmit that frequency 


or to use notch filters to eliminate that frequency and the problem associated with it. 


E. REDUCING THE COMPUTING COST 

Concerning the computing cost of the algorithm, if two full matrix multipli- 
cations are used every time to solve for the spectrum, computing cost becomes a 
factor as each matrix multiplication requires 0(N7) multiplications and additions. 
About half of each inverse matrix is zeros and the solution could be hard-wired more 


efficiently to reduce flop counts. 
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F. GENERALIZING THE SOLUTION 


Having presented the solution method by example, it is necessary to generalize 


the result for any two sampling frequencies. 


Theorem 8.1 Given fs, = Ni where Nj is an even number and fs = No = Ni +1, 


it is possible to form two No x Ng nonsingular matrices A, and A; which solve the 


system of equations A,X =a and AiY = b given by: 


0 1 0 0 00 0 0 0 0 1 
0 0 1 0 0 0 0 0 0 1 QO 
0 0 0 1 0 0 0 0 1 0 0 
0 0 0 0 0 0 0 0 0 0 Q 
00 0 0 1 0 1 0 0 0 Q 
0 0 0 0 0 1 0 0 0 0 Q 
00 0 0 O 1 1 9 0 0 0 
0 0 0 QO Z- OF OQ 0 0 0 
0 0 0 1 0 0 0 90 0 1 0 
0 0 1 0 0 0 0 0 0 0 1 
0 1 0 0 0 0 0 0 0 0 0 
1 0 Q Q 0 0 0 O 0 0 0 
0 1 0 0 0 0 0 0 0 0 
0 0 1 0 0 0 0 0 0 —-1 
0 0 0 1 0 0 0 0 ~1 0 
0 0 0 0 0 0 0 0 0 0 
0 0 0 QO 1 0 -1 0 0 0 
00 0 0 0 1 0 0 0 0 
0 0 0 O 01 -1 0 0 0 
0 0 0 QO 1 0 0 -1 0 0 
0 0 0 1 0 0 0 0 0 -1 
0 0 1 O 0 0 0 0 0 0 
0 1 0 0 0 0 0 0 0 0 
1 0 0 0 0 0 0 0 0 0 
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b2,No 
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The inverses of A, and A; exist and the linear system may be solved X = A, la 


and Y = Aj'b given by: 


Xo 0 0 0 0 
Al 1 1 1 1 
XQ 0 1 1 1 
X3 0 0 1 1 
XN; /2—1 0 0 oO 0 
Xn, /2 = 0 0 0 0 
XN, /2+1 0 0 0 0 
Ny /2+42 0 0 0 0 
XN;-3 0 0 0 +1 
Xn,-2 0 0 —-1 -—I] 
XN,-1 0 —] —] —] 
XN, —-1 -l -—1l -i 
Yo 0 0 0 0 
“4 i 2 a4 
Yo 0: a i 
¥3 0 01 1 
Yn. /2-1 0 0 0 0 
Yn, /2 _ 0 0 0 0 
Yn, /2+1 0 0 0 0 
YN, /242 0 0 0 0 
Yn,-3 00 0 1 
Yn, _92 0 0 1 1 
Yn,-1 Ot 1 A 
Yn, tie. 
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Q1,Ny 
Q2,N, 
Q3,N4 
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an; /2,N1 
AN, /2,N2 

Qn, /2—1,Ne2 


Q3,N5 
22,No 
G1,No2 
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b1N; 
b2 Ny 
b3.N, 
b4,N, 


bn, /2-1,N1 


ON, /2,N; 
ON, /2,N2 


ON, /2-1,No 


63 No 
b2,No 
b1 No 
bo, No 


All methods for solving undersampled signals so far have required at least two 


sampling frequencies. The additional complexity is a trade-off for the additional infor- 


mation gained. In the next chapter, however, a second algorithm is introduced which 


requires only one sampling frequency and simplifies the overall solution considerably. 
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IX. RESOLVING MULTIPLE FREQUENCIES: SECOND 
ALGORITHM 


A. BACKGROUND 

So far the solution method has required that an undersampled signal be sam- 
pled more than once with a different sampling frequency each time. It is, however, 
possible to undersample a signal more than once with a single sampling frequency 
yet still recover the original spectrum. In fact, the groundwork for such a system has 
already been laid in Chapter VIII. Recalling the discussion in Chapter VIII, Section 
C, the problems encountered with the pivotal frequency required that the signal be 


band-limited between 0 and “3: to recover the actual value of bin cs, fet’ It appeared 
2 18 





to be a waste of time but a necessary evil to recover the correct bin information for 
this pivotal bin. Options discussed to remove this apparent constraint were to not 


transmit the frequency or use a notch filter to eliminate it before processing. 


B. USING THE BANDLIMITED SPECTRUM 

Consider a signal s(t) which contains frequency components from 0 to f Hz 
of unknown magnitude and phase. Previously, the solution would have required that 
s(t) be sampled at both fs, = f Hz and fs. = f +1 Hz in order to recover the original 
spectrum. As discussed, the solution would also require that s(t) be additionally 
bandlimited between 0 and f/2 Hz in order to recover the correct value for c 5h. far’ 
Discarding all information gained from the band-limited spectrum, (82) and (83) 
would then be used to solve for the original spectrum. In fact, the solution method 
wastes time because the band-limited spectrum has no aliased components and the 


bin values correspond to the actual values for Zo, Z1,..., Zz. 
: 2 


AY 





1. Recovering the First Half of the Spectrum 


Consider a signal, s(t) containing frequency components from 0 to 10 Hz. If 
the signal is first bandlimited between 0 and 5 Hz and then sampled at 11 Hz, the 


actual frequency coefficients are given by: 


Zo Coll, 
Z\ Chity 
42} _ |} Catz (84) 
23 C311, 
LA C411; 
Z5 C511, 


where Z, = Xn+7Y, represents the actual real and imaginary component at the 


frequency of n Hz and cy ¢,, = Gn,f,, +jbn,f,, corresponds to the value of the nz, bin 


at a sampling frequency of f,;. The additional ZL subscript denotes that this is the 
band-limited version of the signal. 

2. Recovering the Other Half of the Spectrum 

Suppose that on another channel the signal s(t) is not band-limited. If s(¢) is 
sampled at 11 Hz and the DFT is applied and normalized the real components will 


go into bins as: 


40,11 XQ 
ay) X1+ X10 
a211 X2+ XQ 

= 85 
A311 X3 + Xg (85) 
Q411 X4+ X7 
a5 11 As + Ag 


while the imaginary components go into bins as: 


bo 0 
b141 Y= <Yi0 
62.11 Yo-Yo 

= 86 
63.41 ) Cae o ce) 
b411 Y4— Y7 
65.11 Y5 — Y6 


48 








Substituting the values for the Z,’s from (84), (85) and (86) may be re-written 


as: 
X10 1,11 — 41,11, 
XQ A211 — 42,11, 
Xg | = | @311 — 4311; (87) 
X7 A441 — @4,11; 
X6 25 a amen? ay 
and 
Yi0 bri = 0111; 
Yo b211 — 59,11, 
Yg | =— | 5311 — 93,11, (88) 
Y7 baii — 04,11; 
Y6 b5.11 — 65,11; 


Combining (84), (87), and (88) yields: 


ZO COT; 

Zy C111, 

22 C241; 

23 C3.Aly 

Z4 C4l1y 

45, |) >= C511, (89) 
Z6 C5,11 — C5,11z, 
L7 C441 = CA Ai; 
Z8 conj | ¢3,11 — 3,11; 
Z9 C241 = Co 11; 
Z10 Chir Cite 


As can be seen, the solution is much simpler requiring only one sampling frequency 


and far fewer computations. 


C. GENERAL SOLUTION 

Figure 9.1 displays a possible architecture to resolve an undersampled signal 
using the method described. The incoming signal is split into two branches. The 
top branch is the band-limited branch and solves for the first half of the spectrum. 
This information is then combined with the information from the bottom branch to 
recover the undersampled spectrum. 


The solution (89) can be generalized: 
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Theorem 9.1 Let the sampling frequency of a system be f; where f, 1s any odd num- 


ber. A spectrum containing frequencies between 0 and f,—1 Hz may be unambiguously 


recovered provided that a copy of the signal is also band-limited to recover the spectrum 


between 0 and ist This information is then used to recover the other half of the 


spectrum as in: 


CO,fsL 
Clie 
C2. fst 


Cis=1_9, far 
Cis=1_1, far 
Chez) fr 
Cisal fp, ~ Cl fey 
Clertiy pf, Chet-1 fer 
C fet 2 fF, - Cfeal_9 faz 
con) 
C3,fs = C3, fst 
7 _ C2.f sz 
C1 fs as CL yar 


where Zn = Xnt+ Jn and Cam = Onm t+ jonm. 


It would be nice to know that (9.1) puts the Nyquist criterion to rest once and 


for all but a fair treatment requires that the problems with the methods described 


be explored. The next chapter discusses the current limitations of the algorithms as 


well as potentially resolving the same limitations. 
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Figure 9.1: Channel architecture to resolve an undersampled signal utiliz- 
ing a single sampling frequency 
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X. PROBLEMS WITH THE SOLUTION 


A. SIGNALS THAT CAUSE DFT LEAKAGE 

So far only frequencies which directly resolve into single bins without leakage 
have been considered. There is a problem with the solution when the incoming 
frequency components lie in between DFT bins and leakage occurs. When the DFT 
is applied to a sequence with frequencies that correspond directly to DFT bins, the 
result is that the components resolve all of their energy into the bins of the DFT. 


This is not so with frequencies that “straddle” bins. 


B. DFT LEAKAGE EXAMPLES 


Consider a signal containing just one frequency at 2.3153 Hz, given by: 
s(t) = cos(2m(2.3153)t) (91) 


Note that the signal has zero phase. Unfortunately, instead of resolving into just bin 
2 when sampled, the side-lobes of the DFT distribute a portion of the energy into 
many bins. 

1. Algorithm of Theorem 8.1 

Using the algorithm of Theorem 8.1, the signal is sampled at both 10 and 11 
Hz, the DFT and normalized with results shown in Figure 10.1. Displayed are the 
real and imaginary components of the DFT. These DFT results along with (82) and 
(83) yield the spectrum shown in Figure 10.2a for the real part and Figure 10.2b for 
the imaginary part. The result is clearly inaccurate. 

2. Algorithm of Theorem 9.1 

The algorithm of Theorem 9.1 is not nearly as inaccurate for a single frequency. 
Figure 10.1c and Figure 10.1d correspond to both the band-limited and normal spec- 


trums which would be used to solve for the original spectrum. The answer would be 
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exact. The algorithm’s strength lies in the fact that half of the solution is determined 
exactly in the top branch of Figure 9.1. The solution is not totally immune to the 
effects of DFT leakage however. As an example, assume a signal s(t) has frequency 
components in 5 Hz increments from 0 to 10 Hz (0, .5, 1, 1.5, ..., 9.5, 10). The 
sampling frequency is 11 Hz. Figures 10.3a and 10.3b display the real and imaginary 
spectra of the band-limited version of the signal after sampling and DFT’ normaliza- 
tion while figures 10.3c and 10.3d display the non-bandlimited spectra after sampling 

and DFT normalization. Utilizing (90), the original spectrum is calculated and plot- 
| ted in Figures 10.4a and 10.4b while Figures 10.4c and 10.4d show the actual values. 


The solution is inaccurate. 


C. LEAKAGE AND SOLUTION INSTABILITY 

It is not too difficult to discern why the solution methods are unstable. One 
might think that perhaps the leakage from all the other frequencies would somehow 
combine in such a way that the relative magnitudes of the resulting spectrum would 
remain intact. Unfortunately, the symmetrical number system and the linear solutions 
presented both depend on nearly all energy being within the bin that it is expected 
to be in. As an example consider if an incoming frequency at 2 Hz correctly processes 
totally into bin 2 at 10 Hz but somehow processes into bin 3 at 11 Hz. Whereas the 
correct solution is for all energy to be at f = 2 Hz in the final spectrum, this small 


bin error will result in calculating that all energy belongs to f = 8 Hz. 


D. SOLVING THE PROBLEM? 

What is needed ideally is some sort of transformation that places all energy 
in the expected bins without leakage. For instance, for f,; = 10 it is desired that 
f > 9.5 and 0 < f < 0.5 to resolve into bin 0, 0.5 < f < 1.5 to resolve into bin 1, 


and 8.5 < f < 9.5 to resolve totally into bin 9. Clearly no such transformation 


eae 


o4 








exists which will resolve a spectrum exactly without leakage. There is one direction, 
however, that may contain some promise. 

Consider again our signal, s(t), with a single frequency at 2.3153 Hz. Sample 
again at 10 and 11 Hz but this time sample for 100 seconds instead of just 1 second. 
Applying the DFT and normalizing, the DFT now has a bin resolution of 0.01 seconds 
instead of just 1 second. This means that any incoming frequency at any integer 
multiple of 0.01 Hz will resolve exactly into a single bin while all others will cause 
leakage. The results of the DFTs are shown in Figure 10.5. Figure 10.5 shows the 
effect of increasing the bin resolution which has the effect of containing the energy 
within the area of interest. The sidelobes decay well before they “leak” much into any 
of the other bin boundaries. Still, while the energy is now contained in the region of 
interest, there needs to be a method to recover the overall real and complex portions 
of the signal that belong to each integer bin. In the case studied, the goal is to find the 
energy that resolves into the five bins for f,; = 10 Hz and f,2 = 11 Hz, respectively, 
and use the linear equations presented to solve for the spectrum with a bin resolution 
of 1 Hz. 

For the algorithm of Theorem 9.1 the same benefits apply. In fact, increasing 
the bin resolution may hold even more promise. The algorithm holds up fairly well 
when only a little leakage occurs but falls apart when leakage increases. 

It should be emphasized that 10 and 11 Hz have been used as sampling fre- 
quencies throughout for simplicity of presentation. Preliminary study has shown 
that a separation of 1 between the two sampling frequencies gives the best dynamic 
range with the lowest sampling frequencies. It can be shown easily, for instance, that 
the dynamic range of the algorithm of theorem 8.1 is still 0 to 10 Hz for sampling 


frequencies of 10 and 12 Hz. 
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Another aspect of theorem 8.1 is that, unlike the multiplicative nature of 
the dynamic range with the symmetrical number system, the dynamic range of the 
solution is additive. It depends upon the number of linearly independent equations 
that can be set up based upon the number of bins for each sampling eat used. 
For this reason, this algorithm would be most cost-performance efficient using just 


two sampling frequencies as presented. 
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Figure 10.1: Normalized DFT of sinusoid at f = 2.3153 Hz: a) real com- 
ponent for f,, = 10 Hz, b) imaginary component for f,; = 10 Hz, c) real 
component for f,. = 11 Hz, and d) imaginary component for fz. = 11 Hz. 
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Figure 10.2: Spectrum calculated for single frequency at 2.3153 Hz: a) 
real component and b) imaginary component 
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Figure 10.3: Spectra of a signal with frequency components from 0 to 
10 Hz in 1/2 Hz increments: a) real component of sampled bandlimited 
signal, b) imaginary component of bandlimited signal, c) real component 
of sampled original signal, and d) imaginary component of original signal 


og 












Q 1.27 






1 0.6075 0.5267 
3 : 0.4337 - Y 0.8755 6 9 5965 a 
= 0 e 
pe O 0.298 © -0,3565 Ch “0.03045 
a ~0.578 ©) 
0.8457 
0 1 2 3 4 5 6 7 8 9 10 
Frequency (Hz) , 
2) : 
| 
| 
| 
eT | 
- 0 | 
5 ! 





b) 





= 1.558 
S ) @ O) 
ts 9 
a, 0 w @ ") 
‘a 1 -0.09748 9 1253 -0.2063 -0.2397 rh -0.9942 
oa 
0 1 2 3 4 5 6 7 8 9 10 
Frequency (Hz) 


¢) 





a 
g © 0.8323 
ee G S 
3 0.4916 1046. © 
5 © -1,344 Ay 1 55g 04248 -0.5108 
0 1 2 3 4 5 6 7 g 9 10 
Frequency (Hz) 
d) 


Figure 10.4: Spectra of a signal with frequency components from 0 to 
10 Hz in 1/2 Hz increments: a) calculated real component, b) calculated 
imaginary component, c) actual real component, and d) actual imaginary 
component. 
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Figure 10.5: Normalized DFT with 0.01 Hz resolution of sinusoid at 
f = 2.3153 Hz: a) real component fo fs; = 10 Hz, b) imaginary compo- 
nent for fs, = 10 Hz, c) real component for f,. = 11 Hz, and d) imaginary 
component for f, = 11 Hz. 
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XI. CONCLUDING REMARKS 


The essential contribution of this thesis is detailing algorithms or solution 
methods that may hold the key to recovering an undersampled spectrum using one or 
more sampling frequencies. Using the Symmetrical Number System and the Discrete 
Fourier Transform it is a very simple process to resolve one or two undersampled 
frequencies within the dynamic range of the SNS. Additionally, when only resolving 
one or two frequencies, each pairwise relatively prime sampling frequency multiplica- 
tively extends the dynamic range of the system. When resolving an entire spectrum 
of frequencies the solution becomes more complicated but it is still a straightforward 
process. The number of unique bins for each sampling frequency dictates the number 
of linearly independent equations that can be formed and hence dictates the dynamic 
range. Using at most two sampling frequencies demonstrated it has been shown that 
the spectral width that can be recovered using either method is approximately double 
the width predicted by Nyquist. 

Future efforts will attempt to resolve two important questions. First, is it 
possible to apply the solution to a signal which has frequency components that don’t 
conveniently lie directly at bin frequencies? The answer to this question may lie in 
increasing bin resolution and then using a mathematical method to determine fairly 
accurately the total energy that actually belongs to each integer frequency. ‘The 
second question: What happens when more than two sampling frequencies are used? 
The matrix solution for two sampling frequencies is very simple and the solution 
matrix has nice properties. That may or may not be the case when three or more 
sampling frequencies are used. 

In spite of problems that remain with the methods presented, the thesis is 


significant as it presents a direction for resolving an undersampled spectrum. Provided 
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that the leakage problem is solved, either or both of the methods may hold the key 
to greatly extending the ability to sample and resolve extremely high frequencies at 


relatively low sampling rates. 
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